Prom Discrete Hopping to Continuum Modeling on Vicinal Surfaces with Applications 

to Si(OOl) Electromigration 
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Coarse-grained modeling of dynamics on vicinal surfaces concentrates on the diffusion of adatoms 
on terraces with boundary conditions at sharp steps, as first studied by Burton, Cabrera and Frank 
(BCF). Recent electromigration experiments on vicinal Si surfaces suggest the need for more gen- 
eral boundary conditions in a BCF approach. We study a discrete ID hopping model that takes 
into account asymmetry in the hopping rates in the region around a step and the finite probability 
of incorporation into the solid at the step site. By expanding the continuous concentration field 
in a Taylor series evaluated at discrete sites near the step, we relate the kinetic coefficients and 
permeability rate in general sharp step models to the physically suggestive parameters of the hop- 
ping models. In particular we find that both the kinetic coefficients and permeability rate can be 
negative when diffusion is faster near the step than on terraces. These ideas are used to provide 
an understanding of recent electromigration experiment on Si(001) surfaces where step bunching is 
induced by an electric field directed at various angles to the steps. 
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I. INTRODUCTION 

Steps on vicinal crystal surfaces, created by a miscut 
along a low index plane, have long been of great interest 
in both basic and applied researchi High quality crystals 
can be grown through step-flow — the uniform motion of 
more or less equally-spaced steps — and step bunching 
instabilities can create arrays of wide flat terrace sepa- 
rated by closely bunched steps. Other arrangements of 
steps could serve as templates for nanoscale structures 
and devices. 

Most fundamental studies of the static and dynamic 
properties of vicinal surfaces are based on generaliza- 
tions of the classic theory of Burton, Cabrera, and Frank 
(BCF), developed more than fifty years agoi^ This the- 
ory describes the diffusion of adatoms on terraces with 
boundary conditions at steps, which are treated as sharp 
line boundaries. Originally BCF assumed that the steps 
acted as perfect sinks and sources of adatoms so that the 
limiting adatom concentration at the step boundaries al- 
ways reduces to local equilibrium. 

Many extensions and modifications of the BCF theory 
have been suggested to provide a more general frame- 
work for the description of different experiments. One of 
the most important was Chernov's introduction of linear 
kinetic coefficients, which permit deviations from local 
equilibrium at steps^ It was soon recognized that in 
general the kinetic coefficients could be asymmetric^ An- 
other generalization permits step permeability or trans- 
parency, with a term in the boundary condition directly 
connecting the limiting adatom concentration on adja- 
cent terraces^ These generalized BCF models provide a 
mesoscopic or coarse-grained description of surface evo- 
lution with effective boundary conditions at sharp steps, 



and we will generally refer to them as sharp step models. 

Many kinetic instabilities seen in experiments have 
been successfully described from this perspective using 
various combinations of boundary conditions. However 
in general it is not clear how to connect the choices and 
values of the effective parameters in sharp step models 
to the underlying physical processes or how to determine 
the uniqueness of such a mapping. A similar difficulty 
arises in trying to relate "microscopic" parameters in ki- 
netic Monte Carlo simulations of discrete hopping mod- 
els to the effective parameters in a generalized sharp step 
model. Very different microscopic models can sometimes 
seem to give equally plausible mesoscopic descriptions of 
limited sets of experimental data. 

In previous worki we proposed a novel continuum two- 
region diffusion model (CTRM), which gave a rather sim- 
ple and unified description of a variety of current-induced 
instabilities seen experimentally on vicinal Si surfaces. 
We will discuss these experiments in more detail later. 
The model assumes that diffusion rates in a finite region 
around a step could be affected by the different local 
bonding configurations and thus differ from those found 
elsewhere on terraces. By extrapolating the steady state 
concentration profiles to the center of the step region, 
we obtained a mapping of the parameters in the CTRM 
to those of an equivalent classical sharp step model. One 
surprising conclusion was that negative kinetic coefficient 
can arise when the diffusion rate near a step is faster than 
that on the terraces. 

In this paper, we will provide a more systematic way 
of deriving the boundary conditions for the continuum 
sharp step models from a rather general ID discrete hop- 
ping model that permits both asymmetric diffusion in the 
step region as well as step permeability. As discussed by 
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FIG. 1: A schematic plot of the ID potential surface near an 
atomic step. Different D's that have dimensions of diffusion 
constants characterize the hopping rates associated with dif- 
ferent barrier heights. Here, we take the width of the step 
region to be 2a. 

Ghez and Iyer^ such an effective ID model can result 
from averaging over relevant 2D configurations of kink 
and ledge sites on an atomic step. We then use these 
ideas to provide a detailed description of electromigra- 
tion on Si(OOl) surfaces, and find a coherent scenario that 
explains most of the interesting experimental findings. 

II. ID HOPPING MODEL 

The simple ID model that we study is schematically 
shown in Fig. ^ where an atomic step site is surrounded 
by a region of width s with generally different diffusion 
rates, induced by reconstruction or rearrangements of lo- 
cal dangling bonds. As we will see, this difference can 
generate effective kinetic coefficients in a sharp step de- 
scription. 

In general, the width s of the step region with differ- 
ent diffusion barriers should vary for different systems. 
However, we find that the essential physics of the hop- 
ping model is not strongly affected by specific choices of s 
of order of a few lattice spacings a. Thus we analyze the 
algebraically simplest case shown in Fig. ^ w ith half-step 
regions of width a. The more general result, needed in 
the analysis of Si(OOl) experiments, is easily obtained by 
replacing a with s/2, as will become clear later when we 
compare the results of this generic hopping model to our 
previous results. 

We include here two additional physical features of the 
step region as illustrated in Fig. ^ O ne is the possi- 
ble asymmetry in the diffusion processes in the up and 
down half-step regions, described by hopping rates D±/a. 
The D± have dimensions of a diffusion constant, and the 
model is usefully characterized by dimensionless param- 
eters 

R± = D±/D t , (1) 



with D t the diffusion constant on the terraces. 

The other feature we build in is step permeability or 
transparency, characterized in our model by a single pa- 
rameter pk (0 < pk < 1). This can be understood as 
the effective probability in our ID model that an adatom 
hopping to site will encounter a kink site at a given tem- 
perature and thus equilibrate with the solid. This param- 
eter takes account of effects from both kink site density 
and ledge diffusion in a full 2D model. When pk = 1, 
the step site acts as a perfect sink maintained by cither 
enough kink sites or fast ledge diffusion or both, and 
consequently the step site concentration will be pinned 
at equilibrium c eq . In the opposite limit with pk = no 
adatoms are incorporated into the solid. The step site 
behaves like any other terrace site and thus is perfectly 
permeable. We neglect other possible sources of perme- 
ability, including direct hopping over the step region from 
one terrace to another or effects of rapid step motion, 9 
which we believe are less physically relevant for our cases 
of interest. In this section we will also consider only dif- 
fusion fluxes from concentration gradients, and discuss 
the case of driven diffusion from an external field in the 
Appendix. 

We assume that the net flux of adatoms that hop be- 
tween step site and site a can be partitioned into two 
effective contributions: 

J a/i = — [Pk{c eq -6(a)} + (1 -p k ) {6(0) -6(a)}] . 
'a 

(2) 

The first term describes an adatom exchange with prob- 
ability pk involving equilibrated "kink-like" adatoms at 
site with density c eq and the neighboring terrace site. 
The second term involves a similar exchange with prob- 
ability (1 — pk) involving unincorporated "ledge- like" 
adatoms with density c (0) . Only the former involves cre- 
ation/annihilation of adatoms, and the latter is treated 
as a normal diffusion flux that conserves the adatom den- 
sity. 

Similarly, the flux from site —a to is 

J-a/2 = — [Pk{c (-a) - c e J + (1 - Pk) {c (-a) - c (0)}] . 
' a 

(3) 

Since we assume that all the sinks/sources reside only at 
site 0, the net flux of adatoms that hop from site a (—2a) 
to site 2a (—a) takes on the simpler form 

J± 3 a/2 = ±— [6(±a)-6(±2a)]. (4) 
a 

As in many other situations of physical interest, we will 
use the quasi-static approximation to simplify the analy- 
sis. Here we assume that the motion of the step region is 
much slower than the relaxation of the terrace diffusion 
field, so that one can determine the diffusion process on 
terrace sites with fixed positions of the step regions. In 
the quasi-static limit the net change in the number of 
adatoms at each terrace site given by a total flux balance 
must vanish, i.e., dc(x)/dt = for all x = ±a, ±2a ... In 
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particular, at sites ±a, the balance of fluxes is given by 
dc (±a 



= 



dt 



= ± [aJ± a /2 - aJ± 3a / 2 ] 



(5) 



At step site 0, c (0) can be determined by balancing 
the conserved flux terms proportional to (1 — p k ) in Eqs. 
(J2J) and J2J), and is given by 



c(0) 



D+c(a) + D-c(-a) 



D. 



(6) 



III. 



RELATING PARAMETERS IN DISCRETE 
AND CONTINUUM MODELS 



Our task now is to relate the physically suggestive pa- 
rameters R± and p k in the discrete hopping model to the 
kinetic coefficients k± and permeability rate P appear- 
ing in the boundary conditions of a continuum sharp step 
model as in Eq. J7J below. For x > 0, consider a smooth 
continuum concentration profile c (x) that passes through 
the discrete concentrations c(a) and c(2a). (The caret 
distinguishes discrete from continuum functions.) The 
behavior of c (x) at larger x is determined by the phys- 
ical processes on the terraces, but does not need to be 
specified explicitly for our purposes here. 

To make contact with the sharp step model, we rewrite 
the fluxes in Eqs. ©-© in terms of c(x). To that end we 
use a Taylor series expansion to linear order to express 
c (a) = c(a) and c(2a) = c(2a) in terms of c + = c(0 + ), 
the extrapolated limiting concentration as x — > + at the 
sharp step edge in a continuum picture, and its associated 
gradient Vc | + . Similarly, c(— a) and c(—2a) can be 
expressed in terms of c~ and Vc |_, which in general are 
different than c + and Vc + . 

Using Eq. © to eliminate c (0) , and substituting into 
Eq. JSJ, we find that the result can be rewritten in the 
form of a generalized linear kinetics boundary condition 8 
with permeability 



± [AVc |± T^c ± ] = fc± (c± - c eq ) + P (c d 
The kinetic coefficients k± are given by 



k± = 



Di 



Pk 



(R±-l)[l + (l-p k )MY 



where 



M 



R+R— 
(R+ + R-) 



R, 



R. 



R- - 1 R+ - 1 



(7) 



(8) 



(9) 



is symmetric on exchange of + and — . Note in general 
that the ratio of the kinetic coefficients satisfies 



k+ 
k- 



R_ - 1 



R-t 



1 



(10) 



independent of pk ■ The permeability rate P can be writ- 
ten as 



P : 



{\-p k )R + R_ 



Pk {R T -l) (R++R-) 



(11) 



Using Eq. (JSJ in the first factor, we see that P is sym- 
metric on exchange of + and — , and has a finite limit as 
Pk -> 0. 

The final parameter v in Eq. J7J) is zero in our present 
treatment since we used the quasi-static approximation 
to derive Eqs. (0) and ©. In principle, a non- vanishing 
v would arise if we took the flux due to step motion into 
account in the discrete hopping model. However, the 
quasi-static limit is valid in most physical cases of inter- 
est, and thus this additional complication can be avoided. 

Equations 1)711 l|l are the central results in this section. 
As mentioned earlier, we find that the sharp step bound- 
ary condition can indeed be generally expressed using lin- 
ear kinetics with permeability. More importantly, we are 
able to relate the effective parameters in the sharp step 
boundary conditions to the physically suggestive param- 
eters we considered in our generic hopping model. This 
mapping provides a simple way to understand many as- 
pects of clcctromigration phenomena on Si surfaces. 

A notable general feature of these equations is that 
the kinetic coefficients k± are proportional to p k and the 
permeability rate P is proportional to (1 — pk)- The ki- 
netic coefficients characterize adatom exchange involving 
equilibrated solid adatoms at kinks and the adatom gas 
phase, while the permeability rate characterizes adatom 
motion across the step without equilibrating with the 
solid phase. Moreover, the kinetic coefficients k± are in 
general asymmetric on the two sides of the step due to 
the asymmetry of emission and diffusion processes from 
kinks. However, the permeability rate P is symmetric 
since the physical processes of hopping from one side to 
the other without attachment at the step always involves 
the diffusion constants on both sides. 

We now consider some limits of the above general ex- 
pressions to illustrate some interesting features of both 
the kinetic coefficients and the permeability rate. 



A. Impermeable steps, pu — *• 1 

This limit is usually considered in treatments of the 
sharp step model, and we used it to analyze current- 
induced instabilities on Si surfaces^ In this limit the only 
way for the adatoms to go across a step is through attach- 
ment/detachment at kinks, and the permeability rate P 
vanishes. 

The results are conveniently described in terms of the 
attachment/ detachment lengths 



d± = D t /k ± . 
Using Eq. JSJl, these are given by 

d± =a(R ± -l). 



(12) 
(13) 



If we replace the width 2a of the step region in the present 
model by a general value s, we recover exactly the results 
we found earlier— using the CTRM. As shown in the ap- 
pendix, Eq. (|13fl also holds to lowest order even when 
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there is an external driving field that affects processes in 
the step region. 

This shows the general validity of the mapping be- 
tween model parameters that we found earlier by con- 
sidering diffusion driven by a weak electric field. For 
R± > f, corresponding to slower diffusion in the step 
region, the attachment/detachment lengths and kinetic 
coefficients are positive. The kinetics is usually called 
attachment/detachment limited when d± 3> I, with / the 
average terrace width in a uniform step train, or diffusion 
limited when < d± <C /. For R± — 1, d± vanishes, and 
the kinetic coefficients diverge. This forces to equal 
c eq in Eq. J7J) and generates the local equilibrium bound- 
ary condition originally proposed in the BCF model. 

More interestingly, for R± < 1, corresponding to faster 
diffusion in the step region, the attachment/detachment 
lengths and the corresponding kinetic coefficients are neg- 
ative. As we showed earlier^ the sign of the kinetic coef- 
ficients plays a key role in interpreting electromigration 
experiments on Si surfaces, since it determines the stabil- 
ity of a uniform step train for a given current direction. 
This application of our general results will be discussed 
in more detail below. 

In the following, we will characterize the limit pk — > 1 
as defining a perfect sink model, since adatoms can not 
diffuse across a step without attachment /detachment at 
kink sites. As a direct consequence, the two sides of the 
step are decoupled and any change of the microscopic 
rates on one side of the step does not affect the kinetic 
coefficient on the other side. However, as shown above, 
the two sides of the step will in general be coupled for 
Pk < 1 through Eqs. ifTUjl and ijTTJl . and the subsequent 
analysis of step dynamics becomes much more involved. 



d p /a 



B. Very permeable steps, pk 







This limit may be physically relevant at low enough 
temperatures, or slow enough ledge diffusion, or some 
proper combination of both. Here the adatoms 
hop around on the surface without encountering 
sinks/sources in the step region. Thus one expects van- 
ishing kinetic coefficients, but a finite permeability rate, 
and this is indeed what Eq. (JSJ) and Eq. 1|11[) predict in 
this limit. 

As in Eq. (|12[1 . let us define a corresponding perme- 
ability length 



d P = D t /P. 
Then Eqs. QJ and © yield 



2a 



(14) 



(15) 



Similar to d±, the permeability length dp can become 
negative when + R_) < 2, with faster diffusion in 
the step region. Eq. (|15fl is consistent with results 
derived from a continuum phase field models Recently 




FIG. 2: Plot of the dimensionless permeability length dp /a 
as a function of R in the symmetric case for a general pk- 



Pierre-Louis and MetoisiS have argued that negative per- 
meability lengths can explain some novel growth-induced 
instabilities seen during electromigration on Si(lll) sur- 
faces. 



C. Partially permeable steps, < pt < 1 

This is the most general case, where only a finite frac- 
tion of adatoms at the step equilibrate at kinks, presum- 
ably corresponding to intermediate temperatures with 
moderate ledge diffusion. We focus on the simplest sym- 
metric case where D± = D s or R + = i?_ = R in Eqs. 
© _ dJ- The attachment/detachment length becomes 



Pk 



and the permeability length is 



(R-Pk), 



dp = 2a (R - 1) 



(R-Pk) 
{l-Pk)R 



(16) 



(17) 



Equation (|16f) can be understood using the same 
physics as in the perfect sink model. With a finite prob- 
ability pk to encounter a kink, an adatom has to move 
faster in the step region (D s = D t /pk) to maintain local 
equilibrium (d — > or k — » oo) compared with the per- 
fect sink case (D t = D s ). The permeability length in Eq. 
(|17f) is a new feature arising from the possibility that the 
adatoms go directly across the step without equilibrating 
with the solid. This expression shows a fairly complicated 
dependence on microscopic motions characterized by R 
and pk- 

A schematic plot of dp versus R for a given pk is shown 
in Fig. |21 Both d and dp diverge as R — » oo, since all 
motion in the step region vanishes in this limit, dp de- 
creases as R decreases, and stays positive for R > 1. Just 
like the attachment/detachment length, the permeability 
length changes sign from positive to negative as R passes 
through 1, with equal hopping rates in the terrace and 
step regions. However, the permeability length becomes 
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positive again for small enough R when the motion in the 
step region is sufficiently fast (R < Pf.) that the probabil- 
ity of crossing the step without involvement of a kink is 
effectively decreased to a point that it is no longer faster 
than hopping on terraces. 



IV. APPLICATIONS TO ELECTROMIGRATION 
ON SI(OOl) 

We now apply these ideas to current-induced instabil- 
ities on vicinal Si surfaces.— Step bunching is seen on 
Si(lll) surfaces when the electric current is properly di- 
rected normal to the step direction' The uniform step 
train is initially stable when the current flows in the op- 
posite direction. There are three temperature ranges be- 
tween about 850°C and 1300°C where the stable and 
unstable directions are reversed. However, at similar 
temperatures vicinal Si(OOl) surfaces miscut along [110] 
exhibit step bunching from current normal to the steps 
in both directions" Characteristic bunching patterns 
have also been observed for current directed at various 
angles to the steps' 

There is general agreement that in the presence of 
an electric field adatoms acquire an effective charge z*e 
(which includes both electrostatic and a "wind-force" 
contribution arising from scattering of charge carriers), 
and thus experience a field-directed force F = z*eE that 
biases their diffusive motion. However, it is less clear 
what are appropriate boundary conditions in a sharp 
step model for these processes and how they might be 
affected by the electric field and by surface reconstruc- 
tion. The generic hopping model studied in Sec. ITT1 helps 
us shed some light on these issues, and the results can be 
applied to electromigration on both Si(lll) and Si(001) 
surfaces. In this paper we discuss applications to Si(001) 
surfaces. The different Si(lll) instabilities will be dis- 
cussed elsewhere' 

The most notable differences in current-induced step 
bunching on Si(001) and Si(lll) surfaces arise from the 
(2 x 1) surface reconstruction (dimerization) on Si(001), 
which persists up to temperatures of at least 1200°C. 18 
Two characteristic directions on the surface are estab- 
lished by dimerization, either parallel or perpendicular 
to the substrate dimer rows in the orthogonal [110] di- 
rection, denoted by || and _L respectively. Experimental 
evidence suggests that the diffusion along the dimer rows 
is much faster at low temperatures, i.e., D\ 3> -D; 1 ' 

In recent experiments' the bunching behavior was 
studied on dimple geometries, where steps of all orien- 
tations are found. As schematically shown in Fig. 
there are in general two angles needed to describe the 
local geometry of the dimple when the electric field is 
applied, 15 characterized by the angle 9 between direction 
of the electric field and the [110] direction, and the angle 
cp between the field and the local normal to the steps. 

The bunching exhibits interesting angular depen- 
dences. When the current is parallel to the orthogonal 
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FIG. 3: A schematic illustration of the dimple geometry on 
the Si(001) surface, (a) The general view of the dimple with 
the crystallographic directions indicated above. Zooming into 
a given local area of the dimple (the dotted line box) , we show 
the step-terrace configuration with a general direction of the 
electric field, tp is the angle between field direction and the 
local normal to the steps, while 9 is the angle between the 
field and [110]. 9 = 7r/4 corresponds to a field direction along 
[010]. (b) The top view of the dimple when 9 — 0. Zooming 
into the dotted-lined box near the center of the dimple with 
tp = 0, we show a top view of the vicinal surface and a side 
view of the step-terrace configuration. Most of basic physics 
of step pairing and bunching will be illustrated in this simple 
ID geometry with the electric field perpendicular to average 
step position. 



[110] direction (9 = 0), the bunching is observed to be 
strongest in the areas where the current is locally paral- 
lel to the step normal direction {tp = 0), e.g., the dotted 
line box in Fig. [3J3. No bunching occurs in the corre- 
sponding perpendicular directions {tp = 7r/2). However, 
if the current is rotated to 7r/4 off the dimer row di- 
rection {9 = 7r/4), the strongest bunching occurs in the 
areas where the current is perpendicular to the local step 
directions {ip = it/ 2). No bunching is seen in the cor- 
responding perpendicular direction {tp = 0), which in 
the previous case was where the maximum bunching was 
found. In the following discussion, we will first study the 
instabilities for the simplest case as shown in the dot- 
ted line box in Fig. {9 = and <p = 0), and then 
generalize our results to arbitrary 9 and tp. 



A. Domain Conversion and Step Pairing 

Let us begin with the simplest case, where the vicinal 
surface is misoriented in the [110] direction. At equilib- 
rium rather straight Sa steps that run parallel to the 
dimer rows of the upper A terrace alternate with much 
rougher Sb steps that run perpendicular to the dimer 
rows of the upper B terrace' When the field is nor- 
mal to the steps, as illustrated in the boxed region of 
Fig. EJd, the terrace diffusion rates normal to the steps 
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satisfy Df 3> Df. We assume that the dimerization 
persists at least to some extent on both adjacent half- 
step regions around each terrace and will similarly affect 
diffusion rates there. The normal diffusion in the two 
half-step regions around a given step is characterized by 
Df and Df . Taking account of the differences in terrace 
diffusion rates, it seems reasonable to assume at least 
that D? > Df, or 



Di 



Df) (Df-Df)>0. 



(18) 



Special cases of this assumption include classical local 
equilibrium steps where R A = R B = 1 and a symmetric 
step model where Df — Df. 

The assumption here essentially states that the funda- 
mental physics on Si(OOl) surfaces is dominated by the 
alternating reconstruction domains on terraces. Under 
this assumption, it is natural to think of the surface as 
made up of alternating A and B units, where the unit 
a (a = A or B) contains an a terrace together with the 
two neighboring a half step regions. 

We consider here cases where the system is driven away 
from equilibrium only by the electric field. We make 
the usual assumption that permeability does not play 
an essential role in this case and take the limit pk — > 1 
for simplicity. We also neglect evaporation and assume a 
positive effective charge. The perfect sink limit decouples 
the concentration fields on the terraces, and permits a 
simple solution to the steady state diffusion problem in 
terms of exponential functions e^ x , where / = Y-x/ksT. 

In almost all cases of experimental interest, the field is 
sufficiently weak that fs <C fit <C 1, where s and It are 
the width of the step and terrace regions, and we obtain 
piecewise linear profiles for the adatom concentration. It 
is then straightforward to write down the general solution 
for the adatom density in unit a as 



(x) 



+ m\ 



IJ+S 



<X<-- 



2—^—2 



<x< 



t +8 

2 



(19) 



where If is the a terrace width. In the above expression 
the origin is set at the center of the terrace region to take 
maximum advantage of symmetry. It is easy to transform 
the origin to the left atomic step position in accordance 
with the previous discussion on hopping models, and the 
results below will not be altered by any specific choice of 
the coordinate system. 

The mf t can be obtained by requiring continuity of 
concentration and flux at ±Zf/2 and are given by 



m?(l?) 



■ ( </ 



sf(R a -l)/(lf + R a s) 



(If) = c eq lff(l-R a )/(lf + R a s) (20) 



Here R a = Df / Df gives a dimensionless measure of the 
relative diffusion rates in the a unit between the terrace 
and step regions in a direction perpendicular to the step 
direction, and c eq is the average concentration for a uni- 
form step array when / = 0. 

In the quasi-static approximation the step velocities 
are computed by a flux balance. The surface flux nor- 
mal to the step is constant throughout the a unit and is 
exactly given by 



D a r 



l? + . 



eqj 



If + R a S ' 



(21) 



Because of the perfect sink assumption, the fluxes in the 
individual a units on either side of a step are independent 
of each other. Thus the step velocity is easy to compute 
for a given step configuration. 

Consider in particular the initial velocity of step S a in 
a uniform step train (if =lf = l t ). This is given by 



= n (j « - 4) 



= Oc e9 / (l t + s) 



(pf - Df ) l t + (Df - Df ) sR a RP 



(It + R a s) (k + RPs) 
I 



(22) 



where a, (3 = A or B and il is the atomic area. In 
this case the velocities of the two types of steps satisfy 
Vq = —vf. Therefore the initial uniform step array is not 
a steady state. Depending on the direction of the elec- 
tric field, one reconstruction domain expands while the 



other shrinks, creating step pairs separated by the minor 
terrace. With a downhill current one finds double height 
Db steps (consisting of an upper Sb step and a lower Sa 
step with a narrow A terrace trapped in between) sep- 
arated by wide B terraces; the equivalent configuration 
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- Effective step region 

S R S; 



Terrace region 




FIG. 4: A schematic illustration of extrapolation for an effec- 
tive step region. With a step-down current, domain (1x2) 
expands to form an effective terrace region with some typical 
concentration profile c t . On the other hand, domain (2 x 1) 
shrinks to l' and forms an effective step region when combined 
with the two step regions bounding it. ct is extrapolated to 
the dotted-dashed line at x = in the middle of the minor 
terrace which represents the effective "sharp" step. 



with Da steps and narrow B terraces is seen for an up- 
hill current. Experiments show that this field-driven step 
pairing continues until it is balanced at short distances, 
probably by step repulsions, as first suggested by Natori 
et aLf2i in the special case where local equilibrium was 
assumed for all the steps, corresponding to R = R B = 1 
in our model. 



B. Continued Bunching of Paired Steps 

Now let us examine the stability of arrays of such 
paired steps. Assuming that the step pairs (bound- 
aries of the minor domain) with constant spacings persist 
throughout the bunching process, as is shown by exper- 
iments, we can define a symmetric effective two-region 
model that can describe the continued bunching of the 
paired steps. To that end, we treat the minor reconstruc- 
tion terrace together with the two step regions bounding 
it as an effective step region that separates one major ter- 
race from another, as schematically shown in Fig. 0] for 
the case of a step-down current. As shown below (and 
previously discussed^), the bunching behavior is deter- 
mined by the field direction and the sign of the kinetic 
coefficient for the sharp step model associated with the 
effective two-region model defined here. 

In the Appendix we discuss the mapping to a sharp 
step model from a ID hopping model that treats the ef- 
fects of the electric field explicitly, while assuming there 
is a perfect sink at x = in the center of the step region. 
In the present case a minor terrace resides at the cen- 



ter of the effective step region. We can still follow basic 
treatment in the Appendix if we take this into account by 
shifting the origin of the coordinate system by transform- 
ing x — > x— {V + s) /2 and s/2 — * I' /2 + s, where V is the 
width of the minor reconstruction domain. Representing 
the discrete terrace concentrations by a continuum func- 
tion c (x) and Taylor expanding as before about x — + 
— the center of the effective step region — we find that 
the effective sharp interface boundary condition takes a 
form analogous to Eq. I|A5I) : 



D? [Vc|+-/c+] = k a 



-eq ) ) 



(23) 



where Df, a = A or B, is the normal diffusion constant in 
the major terrace. The continuum profile is schematically 
depicted in Fig. 0]for the case of a step-down current. 

Two new features are seen in Eq. (|23(l arising from the 
use of a single effective step region to describe the paired 
steps. First, the major terrace is determined by the cur- 
rent direction in the initial step pairing regime. Second, 
both an effective kinetic coefficient k and an effective 
"equilibrium concentration" c eq appear in the sharp step 
boundary condition. The latter is given by 



l-\f{l' + s) 



(24) 



This can be understood heuristically by noting that the 
effective density in the center c eq should be linearly modi- 
fied by the weak field from its value c eq at the "real" local 
equilibrium step near the lower boundary of the effective 
step region. Similarly, the effective attachment length 
d a associated with the effective kinetic coefficient k a is 
given by 



d" 



£1 

k a 



[R a - l] , 



(25) 



where s — I' + 2s is the width of the effective step region 
and R a = sR a /s is the relative diffusivity in the effective 
two region model defined above. 

Equations i|23fl - i|25|) gi ve the mapping between the ef- 
fective two region model describing paired steps sepa- 
rated by major terraces and an equivalent sharp step 
model. In the steady state where the major terraces all 
have the same width, the surface flux in the sharp step 
model can be obtained from Eq. (|21|l as follows. We re- 
place the parameters c eq , s, and R a by the corresponding 
effective parameters c eqi s, and R a . Clearly I = lt + s rep- 
resents the terrace width in the sharp step model. The 
steady state flux in the sharp step model as a function of 
the terrace width is thus given by 



Jg(l) = D?c eq f 



I 



l + 2d a 



(26) 



Note that a — A or B is determined by the current di 
rection. 

To examine the stability of the above steady state, 
consider a small deviation 8x n = e 



e ut of the n th step 



8 



in the uniform step train, where e n — ee™^ . Here e 
is a small constant and (j> is the phase between neigh- 
boring steps. Then the n th step will move in response 
to the unbalanced flux induced by the changed widths 
of the terraces in front l n = I + £ n {e 1 ^ — 1) and back 
l n -i = I + £ n (l — e~ 1 *). The linear amplification rate 
lu = v n /e n is given by 

Ad a f 

LU = nD?c eq - J —s (l-cos^). (27) 

(I + 2d a ) 

An instability towards step bunching results if d a f > 
with a maximum at <f> = 7T, corresponding to step pairing. 
Note that the direction of the field and the sign of the 
effective kinetic coefficient combine to determine when 
step bunching occurs, as discussed earlier^ 

Using Eqs. H25|) and (|27[) . we sec that to get simulta- 
neous step bunching from current in both directions, as 
seen in experiment, requires 

R A >2+-> R B . (28) 
s 

With a step-down current, the first part of the inequal- 
ity in Eq. I|28|l makes the effective kinetic coefficient for 
the effective step region containing the slower diffusion 
domain positive, which results in a step bunching insta- 
bility. The second inequality in Eq. (|28|) give rise to a 
negative effective kinetic coefficient which produces step 
bunching with a step-up current. Note that this does 
not require negative kinetic coefficients for single steps of 
either kind. 

However, if one assumes the individual steps are at 
local equilibrium, (R A = R B = f), then the kinetic co- 
efficient for the effective step region is negative in both 
cases, and therefore bunching is expected only from a 
step- up current. 



C. Angular Dependence 

It is straightforward to extend the above analysis to a 
general dimple geometry, where the domain conversion 
exhibits interesting angular dependences. A schematic 
view of the experimental dimple is shown in Fig. |3Ji. 
Again, we need to consider the fluxes from the neighbor- 
ing terraces going into the step. Using Eq. J2U, we can 
represent the surface flux as the sum of fluxes along the 
two characteristic directions, 

J/ = cos 9 J s j| + sin 9 J A 1 (29) 

for the front terrace and 

J b = cos 6 J A 1 + sin 6J£ \\ (30) 

for the back terrace of step I in Fig. [3^, where || and 
_L are the directions parallel and perpendicular to dimer 



rows as defined earlier. The angular dependent step ve- 
locity is readily obtained 

4 X) (M = v B cos(29-<p), (31) 

where v B is given by Eq. I|22l) . Eq. (|31|l shows that a 
steady state of paired steps will form on the part of the 
dimple where cos (29 — (p) ^ 0. 

In the following, we will concentrate on two special con- 
figurations that are studied experimentally. 16 The first is 
shown in Fig. |3Jd, where the current is parallel to the 
dimer row direction. In this case 8 = 0, and cosy? char- 
acterizes the angular dependence around the dimple. The 
maximum pairing instability occurs at ip = where the 
current is perpendicular to the step normal direction, and 
no instability in seen at (p — ±n/2. From the previous 
discussion in Sec If V Al and Sec If V fil we can easily see 
that continued step bunching occurs with a maximum at 
ip = 0. 

The other interesting configuration corresponds to an 
upright field parallel to [010] direction in Fig. In this 
case the current is at an angle 9 = 7r/4 from the dimer 
row direction. Hence the angular dependence becomes 
cos(7r/2 — ip) = sin (p. The maximum pairing instabil- 
ity occurs at ip = tt/2, where the current is parallel to 
the steps, and no instability occurs when the current is 
perpendicular to the steps. Again the sharp step model 
corresponding to the steady state can be extracted. The 
subsequent step bunching instability for a parallel cur- 
rent was discussed by Liu et al£2k Their stability analysis 
suggests that step bunching generally occurs for a non- 
vanishing attachment /detachment length d, regardless of 
its sign, when the current is parallel to the average step 
positions. 

The results discussed here are in good agreement with 
experiments. For the angular dependent step pairing, the 
result is consistent with the original analysis by Nielsen 
et al.m& However, our explanation for the subsequent step 
bunching is different. Our analysis provide a simpler sce- 
nario that does not require a tensor character to the ef- 
fective charge. 



V. CONCLUSION 

This paper derives expressions for sharp step boundary 
conditions characterized by linear kinetics rate parame- 
ters k± and P for general BCF type models by appropri- 
ate coarse-graining from a microscopic hopping model. 
k± and P are related to the attachment/detachment ki- 
netics at kinks and to diffusion across the ledges respec- 
tively. In particular, our study shows that both parame- 
ters can be negative when diffusion is faster in the step re- 
gion than on terraces. The possibility of negative kinetic 
coefficients was first suggested by Politi and Villain, 23 
but with no derivation or discussion of any physical con- 
sequences. In the appropriate limit, we recover the map- 
ping previously obtained with the CTRMi Our results 
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also seem consistent with those from phase field models, 9 
while providing a simple and physically suggestive pic- 
ture. 

We then used the perfect sink limit (pk = 1) of the 
general model to analyze current-induced instabilities 
on Si(OOl), where the field represents the only driv- 
ing force away from equilibrium, and found results in 
good agreement with experiment. As we will show else- 
where, this same limit also provides a coherent inter- 
pretation of Si(lll) electromigration experiments, where 
novel step wandering behavior is seeniii Thus we believe 
that the perfect sink model provides a simple and physi- 
cally reasonable description of many electromigration ex- 
periments on vicinal Si surfaces. 
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APPENDIX A: ID PERFECT SINK MODEL 
WITH A CONSTANT ELECTRIC FIELD 

In the generic hopping model discussed earlier, we as- 
sumed that the flux arose only from concentration gra- 



where x is evaluated at discrete lattice sites inside the 
step region. It is easy to write down the solution of Eq. 
l|Al)l as 

c s (x) =c eq + A (e fx - 1) , (A3) 

taking account of the perfect sink at i = 0. Here A 
is a constant that can be determined by continuity of 
fluxes at the boundary between step and terrace region, 
i.e. J s / 2 -a/2 = J s /2+a/2- This gives 

e-f a ' 2 Rct (f + a) - [e-f' a / 2 + eJ a / 2 (R - 1)] c e<? 

~~ g/o/2 t e fs/2 _l\R + e /a/2 _ e -/a/2 ' 

(A4) 

Here dt is the discrete concentration on the terrace site. 

To obtain the sharp interface boundary condition, we 
apply flux continuity J s / 2 + a /2 = J s /2+3a/2, and express 
all the discrete terrace concentrations in terms of the ex- 
trapolated c + and the corresponding gradient Vc | + . In 
the weak field limit that is valid in most experiments, 



dients. We consider here the case where there is a addi- 
tional external driving force from the electric field, and 
take the perfect sink limit pk = 1 used to analyze ex- 
periments. In particular, we examine whether or not 
the kinetic rate parameters in the resulting sharp step 
model could depend on the field as Suga et al. previ- 
ously suggested i2i 

In the absence of the field, the ID potential energy 
surface is similar to that in Fig. ^ where now the site 
at x = is a perfect sink surrounded by a more gen- 
eral region of width s with different diffusion barriers. 
When a weak electric field is applied in the positive x 
direction, the potential energy surface will be modified 
by an amount V = — j Fdx, where F = z*eE, z*e is the 
effective charge. The modification of the potential sur- 
face produces a bias for adatom hopping, which will later 
lead to a convective flux contribution in the continuum 
description. 

The driven flux inside the step region can be written 

as 

J x+a/2 = ^ a ' 2 c s (x) - ^e-f a ' 2 c s (x + a), (Al) 
' a a 

where / = |F| /ksT. The quasi-static approximation 
suggests continuity of fluxes, i.e. J x + a /2 = Jx-a/2, which 
leads to the following equation for the discrete concen- 
tration c s (x), 



(A2) 

I 

we can linearize the exponentials in all of the above ex- 
pressions. To the leading order, we obtain the boundary 
condition as 

±D t [Vc |± T/^] = k± (c ± - c eq ) . (A5) 

where results for both the + and — sides can be given by 
symmetry. Note that the term proportional to / is the 
convective flux induced by the field, which is of the same 
order as the concentration gradient. As mentioned earlier 
in Section IIIII the mapping to the kinetic coefficient is 
independent of the field to lowest order, and is given by 

d± = ^ = ±(R±-l)s, (A6) 

where R± = D t /D s . Equation (| A6|) recovers the results 
we derived earlier from CTRM, and is also consistent 
with the general result in Eq. (15} . 



e~ fa/2 c s {x + a) 



-fa/2 



e^/ 2 c s 



(x 
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